Vortex trapping in suddenly connected Bose-Josephson junctions 
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We investigate the problem of vortex trapping in cyclically coupled Bose-Josephson junctions. 
Starting with N independent BECs we couple the condensates through Josephson links and allow 
the system to reach a stable circulation by adding a dissipative term in our semiclassical equations 
of motion. The central question we address is what is the probability to trap a vortex with winding 
number m. Our numerical simulations reveal that the final distribution of winding numbers is 
narrower than the initial distribution of total phases, indicating an increased probability for no- 
vortex configurations. Further, the nonlinearity of the problem manifests itself in the somewhat 
counter-intuitive result that it is possible to obtain a non-zero circulation starting with zero total 
phase around the loop. The final width of the distribution of winding numbers for N sites scales as 
XN a , where a — 0.47 ± 0.01 and A < 0.67 (value predicted for the initial distribution) indicating 
a shrinking of the final distribution. The actual value of A is found to depend on the strength of 
dissipation. 
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In the past few years, experiments on annular Joseph- 
son tunnel junctions in superconductors [TJ |2] and Bose- 
Einstein condensates [3J 0] have tried to address the 
role of non-adiabaticity in the spontaneous production of 
topological defects, a question that has bearing on early- 
universe cosmology El H]. While a first type of 
experiments [5] have used a temperature quench through 
a second-order phase transition from a normal to a super- 
conducting phase, a second type [3j 0] uses interference 
between initially independent condensates as a mecha- 
nism to trap vortices. In the case of superconductors the 
Kibble-Zurek scaling law [B] relating the probability to 
trap vortices to the quench rate has been tested. Ex- 
periments connecting the independent BECs have simi- 
larly tried to test the role of the merging rate in deter- 
mining the probability for observing vortices in the final 
BEC. Motivated by these experiments we have studied 
numerically the related problem of a ring-shaped Bose- 
Josephson junction array. We would like to stress that, 
while there are similarities between our initial conditions 
and those of the aforementioned experiments, there are 
also qualitative differences that will be discussed later. 
Nevertheless, it is quite conceivable that our findings 
here can be tested in future experiments with ultra-cold 
atomic gases [9]. 

The problem we study here is that of N independent 
Bose-Einstein condensates which upon sudden connec- 
tion become arranged on a ring of weakly coupled con- 
densates. We assume r < £ , where r is the single 
condensate radius and £o is the zero-temperature heal- 
ing length. This condition ensures that no vortices form 
within the individual condensates, leaving us only with 
vortices caused by the phase variation along the ring. At 
t = 0, simultaneous Josephson contacts are made be- 
tween each adjacent pair of condensates. As shown in 
Ref. [10] for the case of two initially independent con- 
densates, a relative phase is quickly established once a 



few condensate atoms have hopped from one side to an- 
other. Each pair of neighboring condensates behaves as 
if a random relative phase <p £ (— tt, tt] is chosen locally. 
However, due to the single- valuedness of the macroscopic 
wave function, there are only N—l independent variables. 
Therefore, in our simulations we choose N — 1 relative 
phases independently, each following a flat distribution 
within the interval (— tt, tt] . The N th relative phase lies 
in the same interval and is determined by the constraint 
that the total phase variation around the ring should 
be lirn (n G Z). From the central limit theorem, we 
know that for N — > oo the distribution of n approaches a 
normal distribution with FWHM = 2.354 criV 1 / 2 , where 
a = 1 / a/12 is the standard deviation for a flat distri- 
bution in the interval (— \, ^]. A key point is to realize 
that the classically stable fixed points correspond to all 
the relative phases being equal (modulo 2tt) to a value 
2irm/N , where to e Z is the winding number or charge 
of the final vortex conguration. To allow our system to 
converge to one of these fixed points we let each link 
follow a semiclassical Josephson equation which includes 
a phcnomcnological dissipation term characterized by a 
single parameter 7. Such dynamics allows the system 
to go through phase slips at individual junctions. Thus, 
generally to ^ n. A number of interesting results are 
obtained: 

(i) The distribution of the final winding number devi- 
ates from the initial distribution for all values of N and 7. 
That final distribution for to is narrower than the initial 
distribution for n, indicating an increased probability for 
low-charge vortex configurations (see Fig. [I]). 

(ii) The width of the final distribution scales with the 
size of the system as AiV Q , where a — 0.47 ± 0.01, inde- 
pendent of 7 and A < 0.67 (normal distribution value), 
indicating a shrinking of the basins of attraction for 
higher winding numbers (see Figs. 2j [3]). For 7 < 3 the 
width of the final distribution shrinks upon decreasing 7 
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FIG. 1: Initial distribution of total phases and final distribu- 
tion of stable winding numbers for N = 10 3 and 7 = 5 for 10 5 




FIG. 2: Red plot shows how the FWHM of the final distribu- 
tion of winding numbers scales with N for 7 = 6. The scaling 
exponent is a = 0.47±0.01 and the prefactor A = 0.55 ±0.05. 
Blue plot shows the scaling of FWHM of the initial distribu- 
tion of total relative phases: a = 0.50 ± 0.01, A = 0.67 ± 0.05 
(color online). 



(see inset of Fig. [3J. 

(iii) If one focuses on initial configurations with n = 0, 
the final distribution of winding numbers in the limit of 
large N is still a Gaussian centered around m = with a 
nonzero spread (see Fig. [4]). This reflects the fact that a 
finite fraction of the initial configurations with zero total 
phase have Josephson coupling energies higher than those 
which correspond to nonzero final winding numbers. 

We start our analysis of the Josephson dynamics 
by stating a theorem: If N BECs with random rela- 
tive phases are coupled by a nearest-neighbor Joseph- 
son coupling on a one-dimensional lattice with periodic 
boundary conditions, a necessary condition to obtain a 
metastable non-zero circulation of winding number 27rm 
is N > 4m, the case of 4m links being marginal. The 
proof is as follows: 

Let us assume that each Josephson junction is de- 
scribed by a two-mode Josephson Hamiltonian: 



H= -E 
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FIG. 3: Prefactor A as a function of 7. Note A < 0.67 for all 
values of 7. Inset shows how FWHM of the final distribution 
of winding numbers scale with 7 for N = 10 3 . 
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FIG. 4: Restricted to configurations J^. <j>i,i+i — 0, this his- 
togram for final winding numbers shows that even in the high 
friction limit one can obtain a non-zero circulation. The above 
simulation uses N = 10 3 and 7 = 50. 



where Ej is the Josephson coupling energy, Ec is the 
charging energy, 4>i^+i is the relative phase between i 
and i + 1 (with i = N + 1 identified to i — 1) and 

rii = Ni — ./V^ ' is the deviation from the equilibrium 

value at condensate i. We assume all 7vf 0) s to be the 
same and initially rt,; = 0, so that ^ i n$ = throughout 
the entire evolution. In the classical limit this Hamil- 
tonian can be mapped into that of coupled rigid pen- 
dula, with the first term denoting the "potential energy" 
and the second term the "kinetic energy" of the pen- 
dula system. Now consider a system with N links and 
a total phase difference of 27rm around the loop. As 
stated earlier, the fixed point corresponding to a circu- 
lation of charge m is given by the configuration where 
all the phases are ip m — 2irm/N (modulo 2ir). Here- 
after, we simplify the notation </>j = (fia+i. To determine 
whether this fixed point is stable we consider a configura- 
tion where (pi — (p m + e; with e, = and £j — > 0. The 
potential energy of this new configuration with respect 
to the fixed point is, up to second order in ej, given by 
AE(ei) = (cosip m ) J2i e i- For the fixed point to be sta- 
ble we should have AE(ei) > 0, which requires N > Am. 
This theorem can equally be applied to a system of XY 
spins coupled by Heisenberg interaction. A corollary is 
that final configurations satisfying N/4 < m < N/2 are 
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unstable. 

For a more generic analysis of the fixed points and 
their basins of attraction we derive from Hamiltonian 
a set of semiclassical equations of motion for the relative 
phases and currents at each junction: 

fa(t) = E c [2ji(t) - j i+1 (t) - j<_i(t))] (2) 

ji(t) = - ski 4>i(t) - jfa(t) (3) 

Here time and energies are expressed in units of EJ 1 
and Ej (H = 1), respectively. It is important to note 
that for cyclically coupled Josephson junctions the vari- 
able canonically conjugate to, say, fa is not (rij — Jii+i) 
but rather the quantity J Q ji(t)dt. We have also added 

a phenomenological dissipative term of the form —^fa 
in the equation of motion for jj while neglecting finite- 
temperature noise. It is important to add this term 
for the system to converge to one of the fixed points. 
From our knowledge of three or more coupled pendula 
we know that the system of equations ^ is chaotic 
[TT] and without any damping would typically explore the 
whole phase space without converging to a fixed point. 
To verify this point, we have investigated the dynamics 
of Lyapunov exponents for the case of N = 3. To en- 
sure that the system is in the Josephson regime we take 
Ec/Ej — 0.01 in all our simulations. We find that 3 out 
of 6 Lyapunov exponents are positive, indicating chaotic 
behavior. We note that the Ohmic nature of the dissipa- 
tive term is only justifed at high temperatures |12j or at 
low temperatures if each condensate lives in a large box 

PI- 

An interesting property of Eq. d2~l is that ^ fa is a 
mathematical constant of motion. However, physically 
the system can still change its winding number by going 
through phase slips at any junction. It will be useful to 
incorporate the above constant of motion by imposing 
the restriction fa S (— tt, tt] only at t — and removing 
it for later times. Of course the physical quantity which 
is observed at the end of the evolution is the Josephson 
current at each junction, which depends on the relative 
phase modulo 2tt. 

In order to generate statistics, we consider a large num- 
ber of different initial configurations, with the relative 
phases and numbers chosen as explained earlier. Equa- 
tions ([2])- (|3j) are then numerically integrated for each 
set of initial conditions. After the average current has 
reached its final equilibrium value, its magnitude equals 
sin(27rm/iV) and the value of the final winding number 
m < N/4 is uniquely extracted. A histogram is then 
plotted for all values of m and its width is recorded. To 
obtain the scaling law we have calculated the width as 
a function of N and fitted it to a function of the form 
\N a . The process is repeated for different values of 7. 

To get a qualitative idea of the dynamics and the role 
of dissipation, we consider a certain class of initial con- 
figurations where fa = ip rn + e while fa — ipi — e/(N — 1) 




FIG. 5: Potential Energy landscape for N = 10 and a certain 
class of configurations: m = 2; fa = 4n/10 + e, cj>j = An/10 — 
e /9; j 7^ *• Winding number zero is the global minimum of 
energy landscape and here occurs at e = 3.67T. 



for 2 < i < N. Given fa, this configuration has the low- 
est potential energy. Fig. [5] shows the potential energy 
for such a configuration as a function of e for N = 10 
and m = 2. The first minimum corresponds to the fixed 
point K2 {fa — <f2 for all i) followed by the fixed point K\ 
(fa — (fx + 2tt ; fa = tp 1 for all i > 1) and so on so forth. 
The global minimum of the energy landscape is the con- 
figuration K with zero winding number. Starting with 
the initial configuration mentioned above, Fig. [5] shows 
the path of steepest descent from K2 to Kq. Starting 
from a local minimum one can characterize the size of 
the basins of attraction by the value e c which e takes at 
the next nearest local maximum. However, one should be 
warned that such an estimate applies only to the specific 
class of initial configurations described above. 

The role played by dissipation can also be elucidated 
by studying that class of configurations. Suppose e> e c i, 
where e c \ is the first critical value of e. The system starts 
at an unstable point and as it rolls down to the fixed point 
with one less winding number, loses kinetic energy due to 
friction. If it arrives at the next stable point with kinetic 
energy less than what is needed to overcome the next 
barrier, then it settles down at the fixed point K m -\. 
However, if it has enough kinetic energy to roll over the 
next barrier then the final winding number would be less 
than (to — 1). A similar role can be envisaged for dissipa- 
tion in the general multidimensional landscape: For large 
7, the system settles down in the nearest valley; for small 
7, the particle may escape the initial basin and lower its 
winding number. Thus low friction enhances, by a mod- 
erate factor, the probability of ending in a low-charge 
configuration, as suggested by Fig. [5] and confirmed by 
Fig. [3] 

For a semi-analytical discussion of the basins of attrac- 
tion we focus on the case of N — 5 (stable m = 0,±1) 
and high friction. Let P(m) be the probability of land- 
ing in a final vortex configuration of charge to, Q(n) 
the initial probability for J^. fa — 2im, and P(m\n) 
the probability to obtain a final charge to conditioned 
to J^ifii — 27m. Below we estimate -P(l) and show 
that P(l) < Q(l). First we note: P(l) = P(1|1)Q(1) + 
P(1|0)Q(0) + P(l| - l)Q(-l). We therefore begin by es- 
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timating P(l|l). The limit of high friction ensures that 
the system follows the path of steepest descent towards 
the nearest stable fixed point. The system always resides 
on the hypersurface S n defined by the constant of mo- 
tion (j>i = 2ixn. Note that, on the surface Si, most of 
the m = 1 configurations correspond to the fixed point 
4>i(t) = 2n/5 (i — 1,...,5), whereas m = can emerge 
from five different fixed points on Si , namely, those of the 
type <j>i(i) = 2ir with <j>j(t) = for all j ^ i (i = 1, 5). 
Likewise, m— —1 is dominated by two sets of fixed points 
on Si: five corresponding to one link having undergone 
a 47r total slip, and ten corresponding to two different 
links each having undergone a 2n slip. Note that, even 
for to = 1 on Si, there are many other configurations 
different from the dominant ones mentioned above e.g. 
4>i = 2ir/5 + 2tt, 4>j = 27r/5 — 27r, and cj>k — 2ir/5 for 
k 7^ i,j (i,j = 1, 5). However, in the limit of large 7, 
those configurations involving many different, mutually 
cancelling phase slips should have negligible probability. 

To calculate the area of the basin of attraction for 
771 = 1 , we define a set of five orthonormal vectors Xi such 
that four of them lie on Si and the fifth vector is perpen- 
dicular to Si . We define our origin on Si by shifting that 
of So along x§ by an amount (fx — 2n/5. The five vec- 
tors are then given by: x% — (l/\/2)(l, — 1, 0,0, 0), £2 = 
(1/^(0,0,1,-1,0), x 3 = (1/720X1,1,1,1,-4), x 4 = 
(1/2)(1, 1, -1, -1,0), x 5 = (1/75)(1, 1, 1, 1,1). 

To obtain the basin boundaries on the four- 
dimensional hypersurface we next write the four inde- 
pendent <j)i's in terms of the in-plane basis vectors Xj's 
(i = 1 , . . . , 4) and transform to spherical coordinates 
(r, 61, &2> $3)- Now, the potential energy is given by 
£ — — Ej J^cos^j and the condition d£/dr — de- 
fines the boundary of the basin of attraction. Shifting 
the origin back to So, the basin boundary for 777 = 1 on 
Si is then given by: 

fx sin (r/i + (fx) + / 2 sm(r/ 2 + ipx) 
+f 3 sin(rf 3 + ipx) + / 4 sin(r/ 4 + ^i) = 0, (4) 



where the various fk — fk(@i, 62, O3) ar e obtained from 
a coordinate transformation. The probability P(l|l) to 
end up with m = 1 having started from any point on 
Si is given by the ratio Ax/ Bx, where Ax is the area en- 
closed by the curve Q on Si and Bx is the total area on 
Si subject to the initial constraints <pi(0) G (— 7r, tt]. Us- 
ing Monte Carlo, we obtain P(l|l) = 0.03. Similarly 
we also calculate P(0|1) and P(0|0) by Monte Carlo, 
both yielding 0.94. Using this second result, the sym- 
metry between 777 = 1 and m = — 1, and the fact that 
P(1|0) + P(0|0) + P(-1|0) = 1, we can also obtain 
P(1|0) = P(-1|0) = 0.03. By contrast, the initial distri- 
butions are Q(0) = 0.6 and Q(l) = Q(-l) = 0.2. Hence 
in the limit of large 7, P(l)/Q(l) = 0.15, which indi- 
cates a shrinking of the initial distribution in favor of 
final zero winding number. Full scale simulations based 
on Eqs. ([2|-([3]) yields for the same ratio 0.14. An exact 
agreement would require consideration of infinitely many 
phase-slip histories. 

In the passing, we would like to note that the above 
analysis holds true strictly in the Josephson regime. Ex- 
periments with fully merging independent BECs [5] or 
the scenario of quasi-condensates in BEC formation as 
envisaged by Zurek [B], always go through an interme- 
diate Josephson regime when adjacent condensates start 
to overlap. However, a complete study of the dynamics 
there would require going beyond the two-mode Joseph- 
son Hamiltonian (JT|) for each junction. This is clearly 
reflected in the outcome of experiments by Scherer et at 
[3] where three independent BECs have been merged to 
form stable vortices in the final BEC. 
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